Cluster Definitions

## Get interactive session ##
#  srun --time=08:00:00 --mem=40G -p int --pty bash

# module purge;source /camp/stp/babs/working/software/modulepath_new_software_tree_2018-08-13;module load pandoc/2.2.3.2-foss-2016b;ml R/3.6.0-foss-2016b-BABS;R;

# sbatch --time=08:00:00 --wrap "module purge;source /camp/stp/babs/working/software/modulepath_new_software_tree_2018-08-13;module load pandoc/2.2.3.2-foss-2016b;ml R/3.6.0-foss-2016b-BABS;Rscript runCM.r" --job-name="rCM" --mem=42G -o rCM.slurm >> commands.txt

# sbatch --time=18:00:00 --wrap "module purge;source /camp/stp/babs/working/software/modulepath_new_software_tree_2018-08-13;module load pandoc/2.2.3.2-foss-2016b;ml R/3.6.0-foss-2016b-BABS;Rscript runCM.r" --job-name="rCM" -p hmem --mem=300G -o rCN.slurm >> commands.txt

Prepare Data

module purge;source /camp/stp/babs/working/software/modulepath_new_software_tree_2018-08-13;module load pandoc/2.2.3.2-foss-2016b;ml R/3.6.0-foss-2016b-BABS;R;

Do Dimensionality Reduction

############################################################################### Perform integrated analysis ##

if (length(Obio@sampleDetailList) > 1) {
    DefaultAssay(OsC_CF) <- "integrated"
} else {
    Obio@parameterList$singleCellClusterString <- gsub("integrated", "RNA", Obio@parameterList$singleCellClusterString)
}


# Run the standard workflow for visualization and clustering This will scale on the most variable features only
OsC_CF <- ScaleData(OsC_CF, verbose = FALSE)


redVec <- names(OsC_CF@reductions)


if (!("pca" %in% redVec)) {
    OsC_CF <- RunPCA(OsC_CF, npcs = Obio@parameterList$singleCellSeuratNpcs4PCA, verbose = FALSE)
}
# t-SNE and Clustering

## Add PCA clusters to data collection ##

if (!("umap" %in% redVec)) {
    OsC_CF <- RunUMAP(OsC_CF, reduction = "pca", dims = 1:20)
}

if (!("tsne" %in% redVec)) {
    OsC_CF <- RunTSNE(OsC_CF, reduction = "pca", dims = 1:20)
}

OsC_CF <- FindNeighbors(OsC_CF, reduction = "pca", dims = 1:20)

Cluster Options

cat(paste(knit(text = chnkVec, quiet = T), collapse = "\n"))

Clusterplot_all_1_0_2

Figure 1A: umap Clusterplot with parameter 0.2 here.

Figure 1B: umap Clusterplot with parameter 0.2 here.

### Cluster_dendrogram__1_0_2

Figure 1A: umap Clusterplot dendrogram with parameter 0.2 here.

Clusterplot_all_2_0_3

Figure 2A: umap Clusterplot with parameter 0.3 here.

Figure 2B: umap Clusterplot with parameter 0.2 here.

### Cluster_dendrogram__2_0_3

Figure 2A: umap Clusterplot dendrogram with parameter 0.3 here.

Clusterplot_all_3_0_5

Figure 3A: umap Clusterplot with parameter 0.5 here.

Figure 3B: umap Clusterplot with parameter 0.2 here.

### Cluster_dendrogram__3_0_5

Figure 3A: umap Clusterplot dendrogram with parameter 0.5 here.

Clusterplot_all_4_0_7

Figure 4A: umap Clusterplot with parameter 0.7 here.

Figure 4B: umap Clusterplot with parameter 0.2 here.

### Cluster_dendrogram__4_0_7

Figure 4A: umap Clusterplot dendrogram with parameter 0.7 here.

Clusterplot_all_5_0_9

Figure 5A: umap Clusterplot with parameter 0.9 here.

Figure 5B: umap Clusterplot with parameter 0.2 here.

### Cluster_dendrogram__5_0_9

Figure 5A: umap Clusterplot dendrogram with parameter 0.9 here.

Clusterplot_all_6_1

Figure 6A: umap Clusterplot with parameter 1 here.

Figure 6B: umap Clusterplot with parameter 0.2 here.

### Cluster_dendrogram__6_1

Figure 6A: umap Clusterplot dendrogram with parameter 1 here.

Clusterplot_all_7_1

Figure 7A: umap Clusterplot with parameter 1 here.

Figure 7B: umap Clusterplot with parameter 0.2 here.

### Cluster_dendrogram__7_1

Figure 7A: umap Clusterplot dendrogram with parameter 1 here.

Clusterplot_all_8_1_3

Figure 8A: umap Clusterplot with parameter 1.3 here.

Figure 8B: umap Clusterplot with parameter 0.2 here.

### Cluster_dendrogram__8_1_3

Figure 8A: umap Clusterplot dendrogram with parameter 1.3 here.

Clusterplot_all_9_1_5

Figure 9A: umap Clusterplot with parameter 1.5 here.

Figure 9B: umap Clusterplot with parameter 0.2 here.

### Cluster_dendrogram__9_1_5

Figure 9A: umap Clusterplot dendrogram with parameter 1.5 here.

Clusterplot_all_10_1_7

Figure 10A: umap Clusterplot with parameter 1.7 here.

Figure 10B: umap Clusterplot with parameter 0.2 here.

### Cluster_dendrogram__10_1_7

Figure 10A: umap Clusterplot dendrogram with parameter 1.7 here.

Clusterplot_all_11_1_9

Figure 11A: umap Clusterplot with parameter 1.9 here.

Figure 11B: umap Clusterplot with parameter 0.2 here.

### Cluster_dendrogram__11_1_9

Figure 11A: umap Clusterplot dendrogram with parameter 1.9 here.

Clusterplot_all_12_2_1

Figure 12A: umap Clusterplot with parameter 2.1 here.

Figure 12B: umap Clusterplot with parameter 0.2 here.

### Cluster_dendrogram__12_2_1

Figure 12A: umap Clusterplot dendrogram with parameter 2.1 here.

Clusterplot_all_13_2_3

Figure 13A: umap Clusterplot with parameter 2.3 here.

Figure 13B: umap Clusterplot with parameter 0.2 here.

### Cluster_dendrogram__13_2_3

Figure 13A: umap Clusterplot dendrogram with parameter 2.3 here.

Clusterplot_all_14_2_5

Figure 14A: umap Clusterplot with parameter 2.5 here.

Figure 14B: umap Clusterplot with parameter 0.2 here.

### Cluster_dendrogram__14_2_5

Figure 14A: umap Clusterplot dendrogram with parameter 2.5 here.

library(AUCell)
plotList <- list()
chnkVec <- as.vector(NULL, mode = "character")


dfHeatmapGenes <- read.delim(
  Obio@parameterList$catRefFile,
  header = T,
  sep = "\t",
  stringsAsFactors = F
  
)

geneSets <- list()

for (i in 1:ncol(dfHeatmapGenes)){
    genes <- as.vector(dfHeatmapGenes[2:nrow(dfHeatmapGenes),i])
    genes <- genes[genes %in% rownames(x = OsC_CF@assays$RNA)]
    geneSets[[names(dfHeatmapGenes)[i]]] <- genes
}


Obio@parameterList[["cat2DplotList"]] <- geneSets



###############################################################################
## Get backdrop

exprMatrix <- as.matrix(OsC_CF@assays$RNA@counts)
#logMat <- log10(exprMatrix+1)

# When using a Seurat object #
logMat <- data.frame(OsC_CF[["RNA"]]@data)

## Load tSNE coordinates ##
cellsTsne <- data.frame(OsC_CF@reductions$umap@cell.embeddings)

## done
FNbase <- paste0("CatScatter_Rankings", VersionPdfExt)
FN <- paste0(Obio@parameterList$reportFigDir, FNbase)
FNrel <- paste0("report_figures/", FNbase)
    

pdf(FN)
    cells_rankings <- AUCell_buildRankings(exprMatrix)

min      1%      5%     10%     50%    100% 

1037.00 1150.45 1648.00 1899.50 3361.00 8212.00

dev.off()

png 2

geneSets <- Obio@parameterList$cat2DplotList

cells_AUC <- AUCell_calcAUC(geneSets, cells_rankings, aucMaxRank=nrow(cells_rankings)*0.05)

## Select thresholds ##


FNbase <- paste0("CatScatterHist", VersionPdfExt)
FN <- paste0(Obio@parameterList$reportFigDir, FNbase)
FNrel <- paste0("report_figures/", FNbase)
            
pdf(FN)
    set.seed(123)
    cells_assignment <- AUCell_exploreThresholds(
        cells_AUC, 
        plotHist=TRUE, 
        nCores=1, 
        assign=TRUE
    )
dev.off()

png 2

## Add data to dfExpr ##

## Plot CatScatters ##
for (i in 1:length(Obio@parameterList$cat2DplotList)){
    HMname <- names(Obio@parameterList$cat2DplotList)[i]
    tag <- gsub("[.]", "_", HMname)
    
    FNbase <- paste0("CatScatterHist_", HMname, VersionPdfExt)
    FN <- paste0(Obio@parameterList$reportFigDir, FNbase)
    FNrel <- paste0("report_figures/", FNbase)
    
    
    selectedThresholds <-  cells_assignment[[i]]$aucThr$thresholds 
    
    if ("minimumDens" %in% rownames(selectedThresholds)) {
        pThr <- selectedThresholds["minimumDens", "threshold"]
    } else if ("Global_k1" %in% rownames(selectedThresholds)){
        pThr <- selectedThresholds["Global_k1", "threshold"]
    } else {
        pThr <- selectedThresholds[1, "threshold"]
    }
    
    if (nrow(cellsTsne) > 15000){
        cex = 0.25
    } else if (nrow(cellsTsne) > 1000){
        cex = 0.5 
    } else {
        cex = 1
    }
    
    
    ## Get AUC matrix ##
    tSNE.df <- data.frame(cellsTsne, cell=rownames(cellsTsne))
    mAUC <- getAUC(cells_AUC)[HMname,rownames(tSNE.df)]
    dfAUC <- data.frame(mAUC)
    dfAUC[["cellID"]] <- row.names(dfAUC)
    dfAUC <- merge(dfAUC, tSNE.df, by.x = "cellID", by.y = "cell")
    
    
    input <- list(
        "x_axis" = "UMAP1",
        "y_axis" = "UMAP2",
        "gene" = HMname
    )
    dotsize <- cex
    
    legendNote <- paste0(
            " The following genes of this dataset are represented in this figure: ",
            paste0(sort(Obio@parameterList$cat2DplotList[[i]]), collapse = ", ")
        )
    
     plotList[[tag]] <- ggplot(data = dfAUC, aes(x=UMAP_1, y=UMAP_2, color = mAUC)
            )+ geom_point( shape=16, size = dotsize
            ) + scale_color_gradient(low="grey", high="darkblue"
            ) + xlab(input$x_axis) + ylab(input$y_axis)  +  theme(
                axis.text.y   = element_text(size=8),
                axis.text.x   = element_text(size=8),
                axis.title.y  = element_text(size=8),
                axis.title.x  = element_text(size=8),
                axis.line = element_line(colour = "black"),
                panel.border = element_rect(colour = "black", fill=NA, size=1),
                plot.title = element_text(hjust = 0.5, size = 12)
            )+ ggtitle(paste0("Category: ", input$gene)
            ) + coord_fixed(ratio = 1) #+ theme(legend.position="none") 
     
     FNbase <- paste0("CatScatter", HMname, VersionPdfExt)
    FN <- paste0(Obio@parameterList$reportFigDir, FNbase)
    FNrel <- paste0("report_figures/", FNbase)
    
    pdf(FN)
        print(plotList[[tag]])
    dev.off()
    ## Create R markdown chunk ##
    figLegend <- paste0(
        "**Figure ", 
        figureCount, 
        ":** Category Scatter showing gene category ", 
        HMname, ". ", legendNote, 
        ". Download a pdf of this figure [here](", FNrel,"). "
    )
            
            
    figureCount <- figureCount + 1 
            
    NewChnk <- paste0(
        "### Category Feature Plot ",HMname,
                "\n```{r CatFeatPlot1_",
                i,", results='asis', echo=F, eval=TRUE, warning=FALSE, fig.cap='",
                figLegend,"'}\n",
                "\n",
                "\n print(plotList[['",tag,"']])",
                "\n cat(  '\n')",
                "\n\n\n```\n"   
            )
          
        
    
    chnkVec <- c(
        chnkVec,
        NewChnk
    )
}


save(OsC_CF,
    file = paste0(
         Obio@parameterList$localWorkDir,
         Obio@parameterList$project_id,
        ".SeuratCF"
     )
)

if (length(plotList) > 3){
    tabVar <- ".tabset .tabset-fade .tabset-dropdown"
} else {
    tabVar <- ".tabset .tabset-fade .tabset-pills"
}

Example Gene categories to help with the clustering decission

cat(paste(knit(text = chnkVec, quiet = T), collapse = "\n"))

Category Feature Plot Enteric_Glia

Figure 15: Category Scatter showing gene category Enteric_Glia. The following genes of this dataset are represented in this figure: S100b, Sox10. Download a pdf of this figure here.

Category Feature Plot Enteric_Neuron

Figure 16: Category Scatter showing gene category Enteric_Neuron. The following genes of this dataset are represented in this figure: Chat, Elavl3, Elavl4, Gad1, Nos1, Ret. Download a pdf of this figure here.

Category Feature Plot Proliferation

Figure 17: Category Scatter showing gene category Proliferation. The following genes of this dataset are represented in this figure: Mki67, Top2a. Download a pdf of this figure here.

Category Feature Plot Apoptosis

Figure 18: Category Scatter showing gene category Apoptosis. The following genes of this dataset are represented in this figure: Acin1, Acvr1c, Aifm3, Akt1, Apc, Bbc3, Bcap31, Birc2, Blcap, Bmx, Bnip1, Capn10, Casp2, Casp3, Casp6, Casp7, Casp8, Cdk5rap3, Cecr2, Cflar, Cideb, Cidec, Clspn, Dedd2, Dffa, Dffb, Dicer1, Dnase1l3, Dnase2a, Dnm1l, Endog, Ern2, Fnta, H1f0, Hmgb2, Htra2, Kpna1, Kpnb1, Madd, Prkcd, Prkcq, Ptk2, Rock1, Satb1, Sharpin, Stk24, Taok1, Tardbp, Top2a, Xkr8. Download a pdf of this figure here.

Documentation

sessionInfo()

## R version 3.6.0 (2019-04-26)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: CentOS Linux 7 (Core)
## 
## Matrix products: default
## BLAS:   /camp/apps/misc/stp/babs/manual/software/r/R-3.6.0-foss-2016b/lib64/R/lib/libRblas.so
## LAPACK: /camp/apps/misc/stp/babs/manual/software/r/R-3.6.0-foss-2016b/lib64/R/lib/libRlapack.so
## 
## locale:
##  [1] LC_CTYPE=en_GB.utf-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_GB.utf-8        LC_COLLATE=en_GB.utf-8    
##  [5] LC_MONETARY=en_GB.utf-8    LC_MESSAGES=en_GB.utf-8   
##  [7] LC_PAPER=en_GB.utf-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_GB.utf-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
##  [1] grid      parallel  stats4    stats     graphics  grDevices utils    
##  [8] datasets  methods   base     
## 
## other attached packages:
##  [1] AUCell_1.6.1                ggtree_1.16.6              
##  [3] sctransform_0.2.1           KernSmooth_2.23-15         
##  [5] fields_10.3                 maps_3.3.0                 
##  [7] spam_2.5-1                  dotCall64_1.0-0            
##  [9] DoubletFinder_2.0.3         scales_1.1.1               
## [11] RColorBrewer_1.1-2          mixtools_1.2.0             
## [13] RMySQL_0.10.20              DBI_1.1.0                  
## [15] DT_0.14                     Seurat_3.1.5               
## [17] knitr_1.29                  forcats_0.5.0              
## [19] stringr_1.4.0               purrr_0.3.4                
## [21] readr_1.3.1                 tidyr_1.1.0                
## [23] tibble_3.0.1                tidyverse_1.3.0            
## [25] ggplot2_3.3.2               dplyr_1.0.0                
## [27] DESeq2_1.24.0               SummarizedExperiment_1.14.1
## [29] DelayedArray_0.10.0         BiocParallel_1.18.1        
## [31] matrixStats_0.56.0          Biobase_2.44.0             
## [33] GenomicRanges_1.36.1        GenomeInfoDb_1.20.0        
## [35] IRanges_2.18.3              S4Vectors_0.22.1           
## [37] BiocGenerics_0.30.0        
## 
## loaded via a namespace (and not attached):
##   [1] R.utils_2.9.2          reticulate_1.16        tidyselect_1.1.0      
##   [4] RSQLite_2.2.0          AnnotationDbi_1.46.1   htmlwidgets_1.5.1     
##   [7] Rtsne_0.15             munsell_0.5.0          codetools_0.2-16      
##  [10] ica_1.0-2              future_1.17.0          withr_2.2.0           
##  [13] colorspace_1.4-1       highr_0.8              rstudioapi_0.11       
##  [16] ROCR_1.0-11            listenv_0.8.0          labeling_0.3          
##  [19] GenomeInfoDbData_1.2.1 bit64_0.9-7            farver_2.0.3          
##  [22] vctrs_0.3.1            treeio_1.8.2           generics_0.0.2        
##  [25] xfun_0.15              R6_2.4.1               rsvd_1.0.3            
##  [28] locfit_1.5-9.4         bitops_1.0-6           assertthat_0.2.1      
##  [31] promises_1.1.1         nnet_7.3-12            gtable_0.3.0          
##  [34] globals_0.12.5         rlang_0.4.6            genefilter_1.66.0     
##  [37] splines_3.6.0          lazyeval_0.2.2         acepack_1.4.1         
##  [40] broom_0.5.6            checkmate_2.0.0        BiocManager_1.30.10   
##  [43] yaml_2.2.1             reshape2_1.4.4         modelr_0.1.8          
##  [46] backports_1.1.8        httpuv_1.5.4           Hmisc_4.4-0           
##  [49] tools_3.6.0            ellipsis_0.3.1         ggridges_0.5.2        
##  [52] Rcpp_1.0.4.6           plyr_1.8.6             base64enc_0.1-3       
##  [55] zlibbioc_1.30.0        RCurl_1.98-1.2         rpart_4.1-15          
##  [58] pbapply_1.4-2          cowplot_1.0.0          zoo_1.8-8             
##  [61] haven_2.3.1            ggrepel_0.8.2          cluster_2.0.8         
##  [64] fs_1.4.2               magrittr_1.5           data.table_1.12.8     
##  [67] RSpectra_0.16-0        lmtest_0.9-37          reprex_0.3.0          
##  [70] RANN_2.6.1             fitdistrplus_1.1-1     hms_0.5.3             
##  [73] patchwork_1.0.1        mime_0.9               evaluate_0.14         
##  [76] xtable_1.8-4           XML_3.99-0.3           jpeg_0.1-8.1          
##  [79] readxl_1.3.1           gridExtra_2.3          compiler_3.6.0        
##  [82] crayon_1.3.4           R.oo_1.23.0            htmltools_0.5.0       
##  [85] segmented_1.2-0        later_1.1.0.1          Formula_1.2-3         
##  [88] geneplotter_1.62.0     lubridate_1.7.9        formatR_1.7           
##  [91] dbplyr_1.4.4           MASS_7.3-51.4          rappdirs_0.3.1        
##  [94] Matrix_1.2-18          cli_2.0.2              R.methodsS3_1.8.0     
##  [97] igraph_1.2.5           pkgconfig_2.0.3        rvcheck_0.1.8         
## [100] foreign_0.8-71         plotly_4.9.2.1         xml2_1.3.2            
## [103] annotate_1.62.0        XVector_0.24.0         rvest_0.3.5           
## [106] digest_0.6.25          RcppAnnoy_0.0.16       tsne_0.1-3            
## [109] graph_1.62.0           rmarkdown_2.3          cellranger_1.1.0      
## [112] leiden_0.3.3           tidytree_0.3.3         htmlTable_2.0.0       
## [115] uwot_0.1.8             GSEABase_1.46.0        kernlab_0.9-29        
## [118] shiny_1.5.0            lifecycle_0.2.0        nlme_3.1-139          
## [121] jsonlite_1.7.0         viridisLite_0.3.0      fansi_0.4.1           
## [124] pillar_1.4.4           lattice_0.20-38        fastmap_1.0.1         
## [127] httr_1.4.1             survival_3.2-3         glue_1.4.1            
## [130] png_0.1-7              bit_1.1-15.2           stringi_1.4.6         
## [133] blob_1.2.1             latticeExtra_0.6-29    memoise_1.1.0         
## [136] irlba_2.3.3            future.apply_1.6.0     ape_5.4